Pairwise Comparisons (GSE137001)

Pathway Heatmaps (Enrichr)

# Get enrichr results and save to file
filename <- file.path(here::here(), "analysis", "avg-expr-GSE137001", "enriched_res.Rdata")

padj_cutoff <- 0.01
lfc_cutoff <- 0.5
dbs <- c("GO_Biological_Process_2021")

if (!file.exists(filename)) {
  enriched_res <- lapply(lfc_dfs, function(df) {
    genes <- df %>%
      filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
      pull(SYMBOL)
    
    return(enrichR::enrichr(genes, dbs))
  })
  save(enriched_res, file = filename)
} else {
  load(filename)
}

Combined Score

# Extract combined scores from each enrichr result
pathways <- c("GO:0009267", "GO:0034198", "GO:0036003", "GO:0036500", "GO:1990928", "GO:0034599", "GO:0034976", "GO:0016236", "GO:0010506")

combined_scores <- lapply(enriched_res, function(res) {
  c_score <- res[["GO_Biological_Process_2021"]] %>%
    filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
    arrange(Term) %>%
    pull(Combined.Score)
  
  return(c_score)
})

combined_scores <- as.data.frame(combined_scores)

row.names(combined_scores) <- enriched_res[[1]][["GO_Biological_Process_2021"]] %>%
  filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
  arrange(Term) %>%
  mutate(Term = str_wrap(Term, width = 60)) %>% 
  pull(Term)

pheatmap::pheatmap(
  combined_scores,
  cluster_cols = FALSE,
  scale = "row",
)

-log10(adjusted p-value)

# Extract combined scores from each enrichr result
pathways <- c("GO:0009267", "GO:0034198", "GO:0036003", "GO:0036500", "GO:1990928", "GO:0034599", "GO:0034976", "GO:0016236", "GO:0010506")

combined_scores <- lapply(enriched_res, function(res) {
  c_score <- res[["GO_Biological_Process_2021"]] %>%
    filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
    arrange(Term) %>%
    mutate(neglog10 = -log10(Adjusted.P.value)) %>% 
    pull(neglog10)
  
  return(c_score)
})

combined_scores <- as.data.frame(combined_scores)

row.names(combined_scores) <- enriched_res[[1]][["GO_Biological_Process_2021"]] %>%
  filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
  arrange(Term) %>%
  mutate(Term = str_wrap(Term, width = 60)) %>% 
  pull(Term)

pheatmap::pheatmap(
  combined_scores,
  cluster_cols = FALSE,
  scale = "row",
)

DEGs and Volcano

state_d2_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
)

state_d4_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
)

state_d6_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
)

state_d9_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
)

state_d12_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
)

state_iPSC_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
)

state_ESC_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
)